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The thermodynamic approach to density functional theory (DFT) is used to derive a versatile 
theoretical framework for the treatment of finite-temperature (and in the limit, zero temperature) 
Bose-Einstein condensates (BECs). The simplest application of this framework, using the over- 
all density of bosons alone, would yield the DFT of Nunes (1999). It is argued that a significant 
improvement in accuracy may be obtained by using additional density fields: the condensate ampli- 
tude and the anomalous density. Thus, two advanced schemes are suggested, one corresponding to 
a generalized two-fluid model of condensate systems, and another scheme which explicitly accounts 
for anomalous density contributions and anomalous effective potentials. The latter reduces to the 
Hartree-Fock-Bogoliubov approach in the limit of weak interactions. For stronger interactions, a 
local density approximation is suggested, but its implementation requires accurate data for the ther- 
modynamic properties of uniform interacting BEC systems, including fictitious perturbed states of 
such systems. Provided that such data becomes available, e.g., from quantum Monte Carlo compu- 
tation, DFT can be used to obtain high-accuracy theoretical results for the equilibrium states of 
BECs of various geometries and external potentials. 

PACS numbers: 03.75.Hh, 67.85.-d, 71.35.Lk 



I. INTRODUCTION 



Our understanding of quantum-degenerate dilute Bosc 
gas systems has increased dramatically since the ex- 
perimental achievement of Bose-Einstein condensation 
(BEC) in ultra-cold dilute alkali gases in 1995 The 
explosive growth of knowledge of these systems has en- 
riched both atomic-molecular-optical physics and many- 
body physics. Landmark developments include the real- 
ization of control over the interaction strength through 
magnetic-field tuning of a Feshbach resonance [3, [f| , the 
generalization to Fermion systems, where the crossover 
between BCS-type and BEC superfluidity has been ob- 
served 0, nonlinear atom optics Q, mixed-phase con- 
densates with different hyperfine states JM , mixtures of 
different bosonic atoms (e.g., Na-Rb) [jfl , mixed Bose- 
Fermi ultra-cold dilute gas systems yd], [l2| , and studies 
where optical potentials have been imposed on conden- 
sate systems, including systems in optical lattices which 
are analogous to condensed-matter systems with cither 
weak or strong correlations (l3| . In the present work, 
advanced methods for evaluating the finite temperature 
equilibrium properties of a single-component dilute sys- 
tem of degenerate bosonic atoms in an external potential 
will be studied. 

A dilute gas of atoms behaves classical ly as long as 
the thermal de Broglie wavelength, Xt = y / 2nh 2 jrak^T , 
is smaller than the mean spacing between atoms, n. -1 ' 3 , 
where T is the temperature, n is the gas density, and 
m is the atomic mass. Quantum degeneracy for a gas 
of bosonic atoms ensues when the atomic wave packets 
begin to overlap, i.e., when At <^ n" 1 / 3 , and a conden- 



sate becomes populated. For a uniform noninteracting 
bosonic system, the critical temperature T c , is given by 
the condition nA^ = C(3/2) ~ 2.612, where £ is the 
Ricmann zeta-function. As the temperature is lowered 
further, the kinetic energy of a uniform noninteracting 
gas decreases, and vanishes in the T — s- limit. For 
a trapped noninteracting gas, all the bosons occupy the 
ground state in this limit, with the zero-point kinetic 
energy density balancing that of the external potential. 

In realistic Bose-condensed systems, interactions can- 
not be ignored, as the interaction energy density is typ- 
ically not small compared with the kinetic-energy den- 
sity (or with that of the external potential). For systems 
which arc dilute enough for the mean spacing between 
atoms to be much larger than the range of the atomic 
potential [3], the interactions can be well modeled with 
a contact interaction with an effective coupling constant 
g = 4irh 2 ao/m, which is proportional to the two-body 
s-wave scattering length, clq (see Ref. [H[). At zero tem- 
perature, the strength of the interactions is quantified by 
the "gas parameter" , nag. In many practical applications 
this parameter is very small, because ao is of nanometer 
scale whereas the density of the alkali atoms at the mo- 
ment of condensation is of order several atoms per cubic 
micrometer. A reasonably good description of such sys- 
tems can be obtained by using the first-order expression 
for the interaction energy. For high accuracy, terms be- 
yond the first order should nevertheless be taken into 
account, because the density at the center of the trap 
increases dramatically as the system is cooled below the 
condensation temperature, and because the relative mag- 



nitude of the leading correction is of order 10-\ 



rather 



than of order nag. 



Furthermore, the value of ao can be 
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made much larger near a Fcshbach resonance, driving 
the system into a strongly interacting regime. At finite 
temperatures, the effects of the interactions are more in- 
volved, as the ratio gn/k^T is also an appropriate mea- 
sure of the strength of the interactions, in addition to the 
gas parameter. If the gas parameter is small, this ratio is 
also small at the transition temperature, but it can be- 
come arbitrarily large upon decreasing the temperature. 

In this paper a finite-temperature density functional 
theory (DFT) approach to treat degenerate bosonic gases 
is developed. It is well known that DFT provides both 
a rigorous conceptual framework and a set of highly- 
accurate practical tools for calculating the ground-state 
properties of interacting electron systems (for an intro- 
duction to DFT see Refs. [H-QJ]). Most calculations of 
the electronic structure of atoms, molecules and solids 
are today carried out using the Hohenberg-Kohn-Sham 
DFT approach introduced in the 1960's [ll Hj. DFT 
has been generalized in many ways, e.g., to treat sys- 
tems at finite-temperature [21( , in time-dependent exter- 
nal fields . (22| - [2^ | superconducting electronic systems J2p 
and systems as diverse as nuclei [26| , classical fluids 
spin density waves [H| and superfluid liquid He 0, 
Another development in DFT is the suggestion of using 
the principles of equilibrium thermodynamics to estab- 
lish the finite-temperature version of DFT as a funda- 
mental thermodynamic representation of the free energy, 
and viewing the ground-state DFT as the T — ► limit 
of this representation [13, Hl|. Here, we follow this ap- 
proach, and apply it to dilute-gas bosonic systems. 

DFT is a method for calculating the energy and den- 
sity distribution of an inhomogeneous system. Within 
the Kohn-Sham approach, it employs a nonintcracting 
reference system which has the same density distribution 
n(r) as the fully interacting system. This noninteract- 
ing system is associated with an effective potential which 
is distinct from the external potential of the interacting 
system. The effects of the interactions may be included 
in a local approximation, which involves a simple inte- 
gral over space, and at each point r, accounts for the 
difference in energy-per-particle between homogeneous 
nonintcracting and interacting systems of density n(r) 
(for electrons, a Hartree term is used to account for the 
long-range part of the Coulomb interactions). Moreover, 
for electrons at zero temperature, this difference in en- 
ergy between homogeneous systems can be described by 
the well-known Wigner interpolation formula [HJ or the 
Gunnarsson-Lundqvist formula (33[, and precise Quan- 
tum Monte Carlo calculations are available [34| . Al- 
though the local density approximation already achieves 
surprisingly high accuracy for many electronic systems, 
the even higher precision required for applications, e.g., 
in chemistry, motivates the ongoing development of more 
sophisticated approaches. Note that DFT is not, in prin- 
ciple, a method for calculating the excitation spectra of 
the systems studied, although the spectrum of the Kohn- 
Sham reference system often fits the spectrum of the 
interacting system quite well (DFT has even become a 



standard tool for evaluating band structure for electrons 
in periodic crystals, although in principle it is only a 
zeroth-order approximation in the context of the meth- 
ods devised for calculating such quantities, such as the 
GW method H). 

In response to the above-mentioned experimental de- 
velopments, several authors developed DFT methods for 
dilute-gas bosonic systems. An early attempt to develop 
a DFT with a high-accuracy Bogoliubov-type treatment 
of the Kohn-Sham system was made in Ref. [36j , which 
employs both the density distribution n(r) and the con- 
densate amplitude $(r) as functional variables. Accord- 
ing to the analysis to be described here, such a high-level 
treatment actually requires the use of three independent 
functional variables, as discussed in Sec. I VII A straight- 
forward application of DFT to boson systems, based on 
the density n(r) alone, was suggested by Nunes [37]] • For 
ground states, i.e., at T = 0, this approach results in 
a modified Gross-Pitacvskii equation, containing terms 
which arc nonlinear in the coupling constant g. It has 
been applied to experimentally relevant regimes 13811 , and 
generalized, e.g., to time-dependent potentials [39|, and 
to address issues particular to strictly one-dimensional 
systems [i(j ■ The possibility of a generalization to finite 
temperatures was noted by Nunes 1371. but it is inferior 
compared to the two-fluid approach (41j (not to be con- 
fused with the Landau two-fluid approach to supcrflu- 
ids), which was already available at the time. Specifically, 
the two-fluid approach achieves improved accuracy by 
treating the condensate component, and the thermal 
component, n— 1$| 2 , separately, with the two components 
subject to different potentials (the interplay between the 
two components can also be studied experimentally [42j ) ■ 
Note, however, that Ref. [4l[ considered large systems 
with very weak inhomogeneities, for which finite-size ef- 
fects are negligible and one may assume local thermody- 
namic equilibrium. For such weakly-inhomogeneous sys- 
tems, application of a sophisticated DFT is superfluous. 
Also note that the two-fluid approach is less accurate 
than the field-theoretic approach in the Popov approxi- 
mation, which was applied at roughly the same time [43| . 

Two different versions of DFT for bosons will be pre- 
sented below, based on the systematic thermodynamic 
approach. One version is based on treating the total den- 
sity, n = and the condensate amplitude, $ = (%/)), 
as two density components. Correspondingly, the Kohn- 
Sham reference system is a nonintcracting boson system 
which has the same n(r) and $(r) distributions, and is 
subject not only to an effective nonintcracting potential 
v n i(r) but also to a fictitious potential ??ni(r) which cou- 
ples directly to the condensate amplitude < I > (r). This ver- 
sion, which may be called $-DFT, reduces to the two- 
fluid approach of Ref. [4l| in the limit of weak interac- 
tions and large systems. It allows for inclusion of ap- 
propriate nonlincar-in-g terms for stronger interactions, 
as well as application to systems with significant inho- 
mogeneities. The second version treats the anomalous 
density A = (ipip) as a third density, resulting in a Kohn- 
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Sham system which is also subject to an anomalous po- 
tential £nj (r) , for which a generalized Bogoliubov-type 
treatment is appropriate. This version is referred to as 
anomalous-DFT or A-DFT, and bears some resemblance 
to the electronic DFT devised for superconducting sys- 
tems [Is)] ■ In the limit of weak interactions, it reproduces 
the Hartree-Fock-Bogoliubov model (a further approxi- 
mation to which yields the Popov model (44j). 

In order to apply a local density approximation for 
the interaction effects, one needs results for uniform sys- 
tems, as discussed above. For dilute Bose gases, some re- 
sults as a function of the density n are available, at both 
vanishing [13, Ell and finite l46M48j temperatures. 
However, application of the advanced DFT versions dis- 
cussed below requires generalization of these results to 
functions not only of the density, but also of the conden- 
sate amplitude for $-DFT and of the anomalous density 
for A-DFT. These generalized interacting systems are 
analogous to spin-polarized uniform electronic systems, 
data for which is in standard use within the local den- 
sity approximation for electrons. The difference is that 
manipulating the condensate amplitude $ or the anoma- 
lous density A requires the use of fictitious potentials, 
whereas the spin density can be modified by subjecting 
the system to a physically realizable magnetic field. As 
the uniform systems which are under consideration are 
fictitious anyway, and the results are obtained by theo- 
retical methods (such as the quantum Monte Carlo work 
mentioned above [34| ) , the realizability or not of the fields 
is of little importance. Obtaining high-accuracy results 
for the generalized uniform systems is beyond the scope 
of the present work, and at this stage we will content our- 
selves with expressions for the interaction effects which 
are valid to first order in g (for A-DFT, this corresponds 
to the Hartree-Fock-Bogoliubov model as noted above), 
with one exception: the leading-order results for A-DFT 
will allow us to deduce next-to-leading-order results for 
<f>-DFT. For dilute gases, the first order approximation to 
A-DFT is wholly sufficient, except for special cases with 
particularly strong interactions, which are realizable near 
Feshbach resonances [1, Hj] . It is also relevant to note that 
homogeneous systems with attractive interactions (nega- 
tive scattering lengths ao) are absolutely unstable at long 
wavelengths, but inhomogeneous attractive systems may 
have metastable dilute-gas states, and have been studied 
experimentally . $-DFT and A-DFT may be applied 
to such systems with the leading-order expressions for 
the interactions, whereas a higher-accuracy local-density 
approach is in principle unworkable, because there can 
be no accurate thermodynamic results for the relevant 
interacting homogeneous system (except at uninterest- 
ingly low densities, where thermal excitations stabilize 
the long wavelength perturbations). The high-accuracy 
methods developed below thus have a range of applicabil- 
ity which is limited primarily to systems with repulsive 
interactions. 

The outline of the paper is as follows. Section II in- 
troduces the general formalism of DFT in the thermody- 



namic language. Section III presents $-DFT: it applies 
the principles of DFT to nonuniform BECs, using the 
total density n(r) and the condensate amplitude $(r) as 
free variables. The presentation includes the Thomas- 
Fermi approximation for the Kohn-Sham system, which 
is applicable when the inhomogeneities are weak, and 
the first order approximation to the interaction energy. 
In Sec. IV, a more general DFT scheme is developed, 
wherein apart from n(r) and $(r), also the anomalous 
density, A(r), is used as a third free variable. In this 
case, the O(g) approximation leads to the Hartree-Fock- 
Bogoliubov system. Here too, the Thomas-Fermi ap- 
proximation is introduced. Section V demonstrates one 
of the advantages of A-DFT, by showing how a result for 
a homogeneous system, which is available with its O(g) 
approximation can be obtained within $-DFT only if 
more complicated higher orders are included. Section VI 
presents a discussion of this comparison, concluding re- 
marks, and suggestions for future research. 



II. FINITE-TEMPERATURE 
DENSITY-FUNCTIONAL THEORY 

The purpose of the present section is to introduce the 
relevant concepts of DFT. The thermodynamic approach 
of Ref. [I?} is followed, and generalized to cases with 
several "density distributions" . This will allow not only 
the total density of particles n(r), but also the condensate 
amplitude $(r) and the anomalous density A(r), to be 
used as free variables, as discussed above. 

The Hamiltonian of the inhomogeneous system may be 
written as H = H nl + Ai?i nt where ff; n t includes all the 
interaction terms, and H nl is a noninteracting (quadratic 
in field-operators) Hamiltonian, for which accurate so- 
lutions are obtainable at an acceptable computational 
cost. A is a continuous parameter specifying the strength 
of the interactions, with A = i for the full interacting 
system, and A = for the noninteracting case. The 
single-particle fields specifying the inhomogeneity, such 
as the potential terms containing the external potential 
v(r) and other fields B(r) are included in H n i, and these 
couple to the densities n(r) and m(r) respectively. In the 
interest of generality, the exact nature of B(r) and m(r) 
will not be specified yet, but as an example one may keep 
in mind electrons in a magnetic field B(r), which couples 
to the spin density (magnetization) m(r) [3~0 ]. 

The focus of the present section is the thermodynamic 
treatment. The statistical-physics problem of obtaining 
the grand potential 17 from the Hamiltonian H will be 
tackled in later sections, where specifics of the Hamilto- 
nian for bosons will be given. In the first subsection here, 
the foundation of DFT will be laid out, by explaining how 
the densities (rather then the external fields) can be re- 
garded as the free functional variables which specify the 
inhomogeneous system. This is a direct generalization of 
Lcgendre transforms, i.e., of the replacement of one free 
variable (e.g., the chemical potential (i) by another (e.g., 
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the total number of particles, N). The second subsection 
explains how, within the DFT framework, the interact- 
ing system can be related to a specific "Kohn-Sham" 
noninteracting reference system, and how the effects of 
interactions can be approximated, based on knowledge of 
homogeneous interacting systems (the local density ap- 
proximation). 




n^N)- 



Legendre Transforms and the Hohenberg— Kohn 
Theorems 



Our starting point uses the grand potential, 
S7([i'(r) — /i, B(r)] , T, A), which depends on the temper- 
ature T and the chemical potential fi, as well as the 
specifics of the Hamiltonian H. The square brackets 
emphasize the functional character of f2, i.e., the fact 
that its value depends on the potential which is itself a 
function of position. The notation also makes explicit 
the fact that the grand potential depends only on 
the difference v(r) — /i, not the values of v(r) and fi 
separately. The derivatives of the grand potential with 
respect to its functional variables are 



n(r) = 



sn 

6v(r) 



m(r) 



sn 

'SB(r) 



(1) 



At this point, these equalities merely introduce notation 
for the derivatives; the fact that n(r) really is the density 
will become evident in the statistical-physics discussion 
of the next section. The different signs used here are a 
matter of convention, and are related to the fact that 
the potential v(r) repels the density n(r), whereas the 
"magnetic field" B(r) attracts the "magnetic moment 
density" m(r), as the magnetic energy density is given 
by — m(r) • B(r). 

We will use the fact that the grand potential fi is con- 
cave in its functional variables, i.e., that when it is eval- 
uated at any two points, [i>i,Bi] and [^2,62] (at fixed 
T > 0, A and /i), and at their midpoint [fi/2, B1/2] with 

t'l/2 = k( Vl + W2 )> B l/2 = 5( Bl + B2 )' tnC mean 01 

the values obtained at the two arbitrary points is strictly 
smaller than the value at the midpoint, |(ili + O2) < 
ili/2 (each one of these "points" is of course a set of 
functions of position, as appropriate for a functional). 
This property, along with others which we shall tacitly 
assume (e.g., differentiability of £1 for finite systems at fi- 
nite temperatures) can be proven by statistical mechanics 
methods (incidentally, Q is also concave in A and T, but 
this will not be used here). The concavity of f2 guaran- 
tees that there exists a one-to-one relationship between 
the potentials and the densities. This corresponds to the 
first Hohenberg-Kohn theorem [TtlUlj]. Thus, a particu- 
lar inhomogeneous system can be identified by its densi- 
ties, n(r) and m(r), instead of specifying the fields v(r) 
and B(r). 

It is convenient to demonstrate this graphically, to- 
gether with the Legendre transforms to be introduced 



FIG. 1: (a) Legendre transform to obtain F(N) from fi(/x) 
(see text), (b) Legendre transform back from F(N) to Q(fi). 
The minimization procedure for the inverse Legendre trans- 
form is demonstrated in the second panel. 



next, using one of the scalar variables of f2, the chemical 
potential /i. The corresponding partial derivative of f2 is 



N 



dji 



(2) 



where N = J drn(r) is the total particle number. The 
one-to-one character of the relationship between \x and 
N follows from the monotonic dependence of the deriva- 
tive of N with respect to the variable fi, see Fig. [TJ Con- 
sider next the combination 0(/x) + fiN, and maximize it 
over all values of /i for a given N. The maximum is clearly 
unique, because the combination is concave in fj,. By con- 
sidering the derivative, one finds that at the maximum (i 
obeys the condition of Eq. This maximum value of 
the combination is called the Hclmholtz free energy, 



F(N) 



max f2(/i) + /iN 



il(n)+(iN\ 



M :iV( M )=JV 



(3) 



The last equality refers to the equivalent procedure of 
choosing \x according to the condition of Eq. ([2]), rather 
than maximizing. The Legendre transform from Q (fj.) to 
F(N) has a simple geometric interpretation (see Fig. [TJ): 
the graph of fl as a function of [i has tangents of slope 
— N, and for a point fl(^t) on the graph, the intercept of 
the tangent line with the vertical axis occurs at F(N) = 

n + fj,N. 

The function F(N), describing the family of tangents 
(intercept as a function of slope) to the curve il(/x), con- 
tains the same information regarding the physical system 
as the original function but is in certain applica- 

tions more convenient. It follows from Eq. @ that the 
derivatives of F are 



dF 



= M 



dF 
dA 



N 



an 

dA 



(4) 



where the last equality represents a derivative with re- 
spect to a variable not involved in the Legendre trans- 
form. As \x increases with N, the function F(N) is con- 
vex (it is concave relative to the other variables, A and 
T). An inverse Legendre transform may thus be applied, 
e.g., by defining 0„(A/') = F(N)— fj,N and identifying the 
grand potential as Cl(fi) = min/v f2 /J (iV) (see the right 
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panel in the figure). The inverse transform differs from 
the original Legendre transform only in signs. 

In the case of functional variables, a geometric inter- 
pretation requires a multitude of "horizontal axes" (one 
for each spatial point) with high-dimensional tangent hy- 
perplanes instead of tangent lines, but the principle is the 
same. The Hohenberg-Kohn free energy of DFT is thus 
introduced through a functional Legendre transform: 

F HK ([n,m],T,A) = Sl([v- l i,B},T,A) 

— J dr {(v — /i) n — B • m} . (5) 

Here the functional variables [v — /i, B] on the right- 
hand-side (RHS) arc determined by maximization, or 
equivalently by requiring the physical condition of Eq. ([TJ) 
(the argument r of the functions is omitted for brevity). 
The functional derivatives of Fhk arc 



A 4 — v ( r ) 



6F, 



HK 



8n(r) 



B(r) 



6F, 



HK 



<5m(r) 



(0) 



The Hohenberg-Kohn free energy, FhkCIt 1 ! 111 ] >T,A), is 
the generalization of the Helmholtz free energy to inho- 
mogeneous systems. 

The inverse Legendre transform allows one to obtain 
the grand potential from the Hohenberg-Kohn free en- 
ergy: 

{l([v-(x,B],T,A) = 
F HK ([n, m] ,T, A) + J dr {(v - M ) n - B • m} . (7) 

Here n and m on the RHS are determined cither by 
Eq. ([6]) or equivalently by minimization. The sec- 
ond Hohenberg-Kohn theorem corresponds to the latter 
statement: the RHS of Eq. (O, when evaluated for an 
interacting system (A = 1) at given external potentials 
v = ^oxt(r) and B — B cxt (r), and minimized over the 
density distributions n (r) and m(r), gives the physical 
value of the grand potential fl at the physical density dis- 
tributions. Although we will make no direct use of this 
minimization principle in the following sections, relying 
instead on Eq. (|6|), its importance in providing both a 
physical picture and an avenue for developing numerical 
algorithms is not to be underestimated. 



B. The Kohn— Sham Equations 

The power of DFT stems from the feasibility of find- 
ing accurate and simple approximations for the compli- 
cated many-body interaction effects in the free energy 
-Fhk- Kohn and Sham exploited the fact that the non- 
interacting effects are much simpler to deal with, and 
nevertheless contain the lion's share of the physics of the 
full system. To introduce the Kohn-Sham scheme, the 



first step is to separate the Hohenberg-Kohn free energy 
into two contributions: 

F HK (Km},r,A = l) = F ni ([n,m],T) + F int ([n,m],T) , 

(8) 

where the nonintcracting free energy F n i of the Kohn- 
Sham system is defined as Fhk in the absence of inter- 
actions, i.e., at A = 0, and the term F nt is defined as 
the difference between the full Fhk at A = 1 and F n j, 
i.e., it contains all of the complicated interaction effects. 
In DFT for electrons, it is standard to further separate 
the interaction term into a simple Hartree long-range in- 
teraction term and an exchange-correlation term which 
is usually considerably smaller, and for which approx- 
imations are sought and employed. As will be clarified 
below, for neutral atoms interacting with short-range po- 
tentials, the "Hartree" (direct) and the "exchange" con- 
tributions are of comparable (often equal) magnitudes, 
and therefore wc proceed with lumping the interactions 
into a single term. 

It follows from Eq. ([5]) that each of the derivatives in 
Eq. ([6]) can also be written as a sum of two terms: 



Uext(r) = v Di (r) - Wi n t(r) , 
B ext (r) = B ni (r) -B int (r) , 



(9) 



where we have used subscripts ext and ni to denote the 
potentials corresponding to the [n, m] densities for A = 1 
and for A = respectively, and the interaction potentials 
are defined as 



^int(r) 



SF ir 



Sn(r) ' 



B int (r) 



SF int 
5m(r) ' 



(10) 



with a convention for the signs which is opposite to 
that of Eq. The external potentials and/or fields 

[w cx t,B oxt ] are known a priori in standard applications, 
whereas the noninteracting potentials [w n i,B n i] (tradi- 
tionally called effective potentials) , which are required to 
reproduce without interactions the same density distri- 
butions [n, m] as in the fully interacting system, are not 
initially known and must be found. Eq. (|9]) immediately 
gives 



v n i{r) = Vcxt(r) + Vint(r) , 

B ni (r) = B cxt (r)+B int (r) , 



(11) 



which gives the noninteracting or effective potentials in 
terms of the externally applied fields plus a contribution 
due to interactions. The system of nonintcracting par- 
ticles in these effective potentials serves as the reference 
system for DFT calculations, and is called the Kohn- 
Sham system. 

Eq. (jlip represents a self-consistent requirement which 
lies at the heart of the Kohn-Sham scheme: given an ini- 
tial guess for the density distributions, and a practical 
approximation for the interaction contribution, this rela- 
tion specifies the potentials for the noninteracting refer- 
ence (Kohn-Sham) system. This reference system may 
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then be solved using the known tools for noninteract- 
ing particles (e.g., the single-particle Shrodinger equation 
with the Fermi-Dirac distribution for the occupations of 
the electrons). The new densities may then be used as 
an improved guess, yielding new values for the nonintcr- 
acting potentials, in an iterative fashion. The iterations 
are stopped once self-consistency has been achieved to 
the desired accuracy. 

It remains to specify the approximation for Fi nt to be 
used. We will limit attention here to local density approx- 
imations (LDAs), of the type suggested (for electrons) by 
Kohn and Sham [2(|. Within this approach, the inter- 
action term is approximated by using the properties of 
uniform interacting systems: 



dr/i„ t (n(r),m(r)) 



(12) 



Here /i n t (n, m) is the contribution of interactions to the 
Hohcnberg-Kohn free energy of a uniform system with 
densities n an m, calculated per unit volume. With this 
simple expression for the interaction term, the functional 
derivatives defining the contribution to the potentials, 
Eq. (TTU)) , can easily be taken: 



Vint(r) = ^^(n(r),m(r)) , 
on 

B int (r) = -^^(n(r),m(r)) 
am 



(13) 



Knowledge of fi n t(n, m) comes from outside of DFT. The 
uniform system is much simpler than the non-uniform 
system in principle, but evaluation of the many-body ef- 
fects even in the uniform case can require sophisticated 
techniques. For example, for electron systems, quantum 
Monte Carlo techniques have been employed, as already 
noted. Once the results are found, the function fi nt (n, m) 
can be tabulated or otherwise efficiently represented. The 
results of the sophisticated calculations for uniform sys- 
tems are thus imported, using DFT, as input for the cal- 
culations of inhomogeneous systems. 

It is of interest to note that the thermodynamic deriva- 
tion used here is constructive. For example, it imme- 
diately gives the exact relation i*} nt = J Q dA (dF/dA), 
with the integrand (dF/dA) equal to d£l/dA — (Hint), 
which in the context of electrons has been called the adi- 
abatic connection formula (50j , and has been derived via 
a much less direct route. For weakly interacting bosons, 
it is appropriate to approximate the integrand here by 
its noninteracting value at A = 0. As we will see below, 
this yields a particularly simple approximation for Fi n t, 
which is again local, i.e., of the form of Eq. (JTSJ) . 



separate field (the condensate field $), in addition to the 
density, will be developed. A system of identical bosonic 
atoms of mass m in an external potential i> ex t(r) can be 
described, in second-quantized notation, by the Hamilto- 
nian 



H = drft I - 



h 2 V 2 

2m 



+ «ext V' + Hint , 



(14) 



where the interaction term involves a two-body inter- 
action potential V(r — r'). This potential has a hard- 
core repulsive form at small interatomic separations, and 
a long-range attractive van-der-Waals form outside the 
core. As explained in the introduction, for a dilute gas, 
with the typical distance between atoms much larger 
than the range of the potential, only the s-wave scatter- 
ing contribution is significant, and the interaction can be 
fully characterized by a single parameter, eto, the s-wave 
scattering length. 

The scattering length ao is typically of the order 
of nanometers, while the typical distance between the 
atoms, n -1 / 3 , in experiments, is of order hundreds of 
nanometers. One may therefore use a Hamiltonian with 
a point interaction, 



— I drip^^^ptp 



(15) 



where a high-momentum cutoff hk c is implied, i.e., no 
attempt to describe the components of the field ip on 
length-scales as small as the range of the interaction 
potential is made. Note that physical quantities de- 
rived from this Hamiltonian depend on both the inter- 
action strength g and the cutoff k c . For example, the 
s-wave scattering length is related to the parameters in 
the Hamiltonian by ao — gm/Aith 2 only to leading or- 
der in g, with corrections of order k c a^, which will be 
assumed small. This undesirable feature may be avoided 
by using a short-range pseudo-potential |45l [5lj , an op- 
tion which will not be made explicit here, but is neces- 
sary when large values of ao are encountered (Fcshbach 
resonances). Note that different forms of H lnt are legit- 
imate within DFT as developed below, and are associ- 
ated with different interaction contributions F- mt . Thus, 
when the pseudo-potential form of H- m t is used, and the 
corresponding changes are made in F mt , all of the DFT 
expressions to be derived below will remain valid (expres- 
sions for Fi n t beyond the leading order are not included 
in the present work). Furthermore, one may include, e.g., 
three-body interactions, simply by modifying Fiat appro- 
priately. 



III. DENSITY FUNCTIONAL THEORY FOR 
BOSONS — $-DFT 



A. The Grand Potential and the Free Energy 



In this section, a version of DFT adapted to bosonic 
systems, in which the condensate will be treated as a 



For a bosonic system coupled to a particle reservoir 
at chemical potential [i and temperature T, the grand 
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potential may be written as 



B. The Kohn— Sham Equations 



ft ([v — fx, 7], if] , T, A) = —k^T In Tr cxp 



"knT 



(16) 

where the trace is over the full many-body Hilbert space. 
Fictitious potential fields, f?(r) and f?*(r), which break 
the particle-number conservation symmetry, have been 
included here in order to couple to the condensate fields, 
$(r) and 3>*(r) which will be introduced shortly. The 
grand-canonical Hamiltonian, ff W = H — /iN with N 
the number operator, is 



Ai? int . 



h 2 V< 

2m 



+ v — /i)?/> — rj'ijr — 7]*tjj 



(17) 



It will be convenient to treat rj{r) and r]*(r) (and simi- 
larly $(r) and $*(r), see below) as independent, and to 
set them equal to the complex conjugates of each other at 
the end of the calculation. Clearly, the physical fictitious 
fields vanish, %xt(r) = ^oxtl 1 ") = 0> but the noninteract- 
ing or effective fields, rj n i(r) = ryi nt (r), may be significant. 

The statistical-physics definition of Eq. (fT6|) fulfills all 
the thermodynamic requirements assumed in the previ- 
ous section. Specifically, it follows directly from Eq. (jT6]) 
that ft is concave Its functional derivatives are as 

described in Eq. ([1} , where we can now identify the den- 
sity as 



n(r) 



sn 

6v (r) 



= (^(r)V3(r)) , 



and the condensate field 



$(r) = - 



on 

dri*(r) 



(18) 



(19) 



with the corresponding expression for $* implied. The 
entropy and interaction energy are given by 



S 



on 

df 



on 

dA 



(20) 



The principles of DFT detailed in the previous section 
may now be applied, with the fictitious potential and 
the condensate field replacing the "magnetic" terms B (r) 
and m(r). The Hohcnberg-Kohn free energy is thus 

F HK ([n, $, $*] ,T, A) = n ([« - ft] , T, A) 

- J dr {(v- fi)n-ri*<f>-ri$*} , (21) 



and its derivatives are 
SF HK 



5n 
dT 



-(v-fi) 



8Fi 



HK 



= -s 



dA 



6$ 



(22) 
(23) 



We next apply the Kohn-Sham approach, based on the 
partition in Eq. of -Fhk into a term describing a nonin- 
teracting reference system and an interaction term. The 
nonintcracting reference system, i.e., the Kohn-Sham 
system, is described by the grand-canonical Hamiltonian 
(we drop the \i superscript to simplify notation) 



H ni = j dr{ft {- 



h 2 v 2 

2m 



■+v ni - f i)^-r ]ni ft-r ] *^}, (24) 



with the nonintcracting (or effective) potentials given by 
Eq. (|11[) . The field operator ip(r) may be written in terms 
of the condensate field $(r) and a residual operator field 
0(r): 



V>(r) = $(r) +0(r) . 



(25) 



The requirement (4>(r)) = 0, cf. Eq. ([!§)) , is associated 
with a modified Gross-Pitaevskii equation, 



h 2 V 2 

2m 



fj, $ = T] n 



(26) 



This condition leads to the vanishing of the linear-in- 
(and similarly, in <f>) terms in the Hamiltonian, 



H n \ = H t h 

At 



"I - H corv 

h 2 v 2 



2m 



(27) 

and hence to (<f>(r)) = 0. A partial cancellation of the 
term involving the condensate field has occurred here, 
and we have introduced notation separating "thermal" 
and "condensate" parts. 

The Schrodingcr equation associated with the thermal 
part of the Hamiltonian is 



ft 2 V 2 
2m 



(28) 



for the single-particle wave functions ipj and their eigen- 
values £j (i.e., the single-particle energies, measured from 
the chemical potential). Using eigenstate creation and 
annihilation operators, 

4>(r) = ]T & (r) a, , 0+ (r) = ]T <p* (r) fit , (29) 

3 3 

one may rewrite this effective noninteracting many-body 
Hamiltonian as 



(30) 



It then becomes straightforward to evaluate the 
statistical-mechanical properties of this noninteracting 
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Kohn-Sham system. The grand potential, Eq. (jT6|) , also 
separates into two parts, 

^ni (Kii - U, Vni, Vni] , T ) = fi th ([«ni - (A > T ) + 

^COIl 

[v ni - /U,7?m,r?^] , (31) 

where 



^th ([fni - /z] ,T) = fc B T V" In ( 1 - exp ( - 



3 



(32) 

(the requirement that the chemical potential be lower 
than the ground state of the Schrddingcr equation, 
mim, £j > 0, is manifest here), and 

ficonflVni ~ «ni, = - ^ / dr (r/ni$* + 7?*;$) . 

(33) 

The Shrodinger equation, Eq. (f2"5| , determines the eigen- 
values in n t h, and the Kohn-Sham form of the Gross- 
Pitaevskii equation, Eq. (|26[) determines the condensate 
field in O con . Note that the effective noninteracting po- 
tential v n i appears in both, whereas the fictitious poten- 
tial r] n i appears only in the latter, and that f2 C on does not 
depend on the temperature T. 

Turning to the functional derivatives, one finds that 
the density distribution, Eq. (fl"8)) . becomes 



n(r) = n th (r) + n con (r) 



E 



expfo/T)-l 



+ $*(r)*(r) , (34) 



together with <I> = — (5f2 C on/<5?7ni- These relations are not 
only obvious from Eq. (|25[) , but can also be derived from 
Eq. (|26|) . Explicitly, one takes its variation and multi- 
plies by $* to obtain 

+ «ni - /«) 5* + - $*^ni , (35) 

where the first term may be identified as n ni S(fr, and then 
5Q con = Jdr (-$*<fy ni - $577*; + **$tfvni) (36) 
follows. 

The Hohenberg-Kohn free energy is 

F ni ([n, $•] , T) - F th ([n - , T) + F con [$, $*] , 

(37) 

with 



Fa, ( [n th ] , T) = fc B r ^ In ( 1 - exp ( - 



dr n t h («ni - A f ) 



(38) 



and 



F con [*,$*] = ydr$* $ • (39) 



As is generally the case with Legcndre transforms, 
the RHS of Eq. (|3"8")l is evaluated for the potential 
7j n i (r) which corresponds to the given density nth(r), 
and it is difficult to make it more explicit. How- 
ever, use of Eq. ([26]) has yielded a significant sim- 
plification in the condensate term, defined as -F con = 
Ocon- / dr ((v ni - n) n con - r) nl $* - r) ni $), resulting in 
the explicit form of Eq. (f3T)l) . The rule for Legcndre 
transforms of sums such as f£ n i = f2th + Q C on is that 
each term can be transformed separately, Q t h into -Fth 
and fl con into F con , but the sum must be evaluated 
as F nl [n] = F t h [nth] + F con [n con ] with the conditions 
n = n t h + iicon and 5F t h/5n t h = SF con /5n con implied. 
In the present case the contribution of the condensate to 
the density, n con = <£>*<!> is known in terms of the con- 
densate amplitude, which is itself a free variable, and no 
implicit relationship remains to be evaluated. In other 
words, the fact that 5tt con /5v n i = <f>* ( f> is trivially re- 
lated to (5f2 C on/<5'7ni = — $ plays a significant simplifying 
role, resulting in F con depending only on $ and $*, and 
F t h depending only on n — 

In summary, the Kohn-Sham equations for a system 
of bosons are Eq. (|26|) for the condensate field, Eqs. I|28[) 
and (J2U) for the density, and v n \ = v cxt + v- m t and rj n i = 
r]i n t for the effective potentials, from Eq. (fTTj) . For an 
LDA, we have 

Vint = dfint/dn , ryint = -df int /d$* , (40) 

where specific expressions for the interaction energy den- 
sity, /i n t(n, $*), will be suggested below. Once these 
Kohn Sham equations have been solved, the grand po- 
tential may be evaluated as 



Q = Cl ai + F int - / dr {{v ni - v cxt ) n - {r) ni <I>* +n ni <f>)} 
= £fc B Tln(l-exp(--g-))+ (41) 



dr { (/ int - nu int ) + - (r? in t$*+Ct $ ) 



The integral here is a generalized subtraction of the 
double counting of the interaction energy included in 
the single-particle energies, as customarily occurs in 
Hartrcc-likc schemes. 

Note that for Fermions there is no condensate term, 
and in the low-temperature limit, F n i is simply the ki- 
netic energy K. The temperature and entropy can be 
thought of as a correction which is necessary at finite 
temperatures. For noninteracting bosons, one still has 
Fth = K + TS, but both terms vanish as the temper- 
ature is lowered, and for large condensates the zero- 
point kinetic energy, F con , may also be negligible. In 
such cases, one has no significant contribution to the 
Hohcnbcrg-Kohn free energy from the Kohn-Sham sys- 
tem, and F ~ F mt in the low-temperature limit. 
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The Thomas-Fermi Approximation 

Many of the relevant Bose-Einstcin condensate sys- 
tems studied experimentally involve a large number of 
bosonic atoms, in the thousands or millions, in a smooth 
external potential. In such cases it is appropriate to in- 
troduce the Thomas-Fermi approximation [ 1 51 ] (adapted 
from many electron systems), which takes the density of 
single-particle states in phase space to be (2wh)~ , and 

2 

uses the classical relationship e — ^— + u n i(r) — \i. The 
local density of states is thus approximated as 



(a) 0.0 



d{e,r) 



d 3 p 

(27Tft) S 



6 (e - «ni(r) + fx) 



2^ - Vni (r) 



(42) 



TOi/2m (e — v ni (r) + jj) 



2tt 2 H 3 

and the overall density of states is given by 



d(e) 



5(e - Ej) 



J drd(e,r) 



(43) 



With this approximation, there is no need to solve the 
Schrodinger equation, Eq. (|28[) . which is the step which 
is most significant in terms of computational resources. 
Expressions such as Eqs. (|32p and (f3~4"|) are then evalu- 
ated as simple integrals over the corresponding density 
of states, Eq. (|43]) or (|42]) respectively. For example, the 
nonintcracting grand potential, from Eqs. (|32[) and 
becomes 



^th = k^TXr, 



drf 



Vni - H 



(44) 



where At is the thermal dc Broglie wavelength mentioned 
in the introduction, and 



f(x) = -= / q 2 dq In 1 - e 



00 -lx 

~ ~ J5J2 



1=1 



(45) 

with x = (v n { — /i)/fceT, and q a scaled momentum vari- 
able. The function f(x) is known as the polylogarithm 
or de Jonquiere 's function, and is plotted in Fig. [5] It 
varies from —£(5/2) to as a; is varied from to 00 (the 
Riemann zeta-function evaluates to £(5/2) = 1.341...). 
Its derivative (also plotted), which varies from £(3/2) = 
2.612 ... to 0, determines the density, 



1 \ \ — 3 tl I 

n th (r) = X T f [-^r- 



Its Legendre transform (cf. the figure again) is 
f(u) = max{/(i) - ux} , 



(46) 



(47) 



with u = AyTith, the dimensionless density, and this de- 
termines the nonintcracting Hohenberg-Kohn free energy 



fth([ntij , T) = k B T\r 3 J dr /(A T n th ) 



(48) 




FIG. 2: Relationships for the thermal components according 
to the Thomas-Fermi approximation: / is the scaled grand 
potential, x is the scaled effective potential, u is the scaled 
density, and / is the scaled Hohenberg-Kohn free energy, (a) 
The function f(x) (full line), and its derivative, u(x) (dotted 
line), (b) The Legendre transform f(u) (full line), and its 
derivative -x(u) (dotted line). Note that the derivatives are 
simply inverse functions. 



The Thomas-Fermi approximation is appropriate for 
systems with gradual inhomogencities. It may be ap- 
plied to the condensate component as well, by simply 
dropping the gradient term in Eq. (|26p . which amounts 
to neglecting the zero-point energy of the condensate, 
F con ~ (this corresponds to f2 CO n — / dr(v n i — //)$*<]> in 
the above notation). As mentioned in the introduction, 
when all the finite-size effects due to the inhomogeneities 
in the system are indeed negligible, it is appropriate to 
use a local-equilibrium approach, with the density n(r) 
at each position taken as that which for an infinite system 
would correspond to the given local value of the chemi- 
cal potential, fi — w ox t(r). Applications of DFT to such 
situations approach the local-equilibrium results. For $- 
DFT, the condensate amplitude $(r) relaxes to the value 
corresponding to an infinite system of density n(r), and 
there is thus no point in including it as a separate func- 
tional variable. 



10 



C. Interaction Effects 

In order to complete the DFT description, an 
approximate description of the interaction term, 
Fint([n, T) must be specified. As noted at the 

end of Sec. II, the simplest approximation is obtained 
by equating the integrand in the adiabatic connection 
formula with its value for noninteracting bosons: 



^ int ([n,$,3>*],T) 



dr - ( 2n 2 



($*$)' 



(49) 



or / int = (g/2) (2n 2 - ($*$) 2 ) in the notation of Eq. 

([T2]) . The factor of 2 in the brackets comes from count- 
ing both the direct and the exchange contributions, and 
the subtraction comes from the fact that exchange is not 
relevant to the condensate's interaction with itself. From 
Eq. g0]), this leads to 



v int = 2gn, r? int = g , 



(50) 



or to 



v ni = v cxt + 2gn , r] ni = g ($*$)$ . (51) 



Using the latter in Eq, (|26p gives the Gross-Pitaevskii 
equation, 



h 2 V 2 

2m 



u ext - n + 2gn th + $ = . (52) 



As mentioned in the introduction, the interaction terms 
here differ from the simple 2gn appearing in the effec- 
tive potential, due to the absence of an exchange contri- 
bution to the condensate-condensate interactions. The 
first-order approximation of Eq. (|49[) . with /; nt quadratic 
in n, $ and <&*, leads to simplification of Eq. (|41 [) for the 
grand potential, resulting in O = O t h — / fintdr, where 
the subtraction of the double counting of the interaction 
energy is explicit. This simplification is not as dramatic 
as it may seem, as the subtraction can only be evalu- 
ated after the Kohn-Sham system of equations has been 
solved (cither with or without the Thomas-Fermi approx- 
imation for the thermal cloud and for the condensate). 

Note that Eq. (j52")l is identical to the G-P equation 
derived from the field-theory approach in the Popov ap- 
proximation 44 1 . 

The present level of description, with the Thomas- 
Fermi approximation for the density of the thermal com- 
ponent, Eq. (|46[) . reproduces the two-fluid description of 
finite-temperature Bose-Einstein condensates [4l| men- 
tioned in the introduction. The present $-DFT provides 
a route for improvements in this description, based on 
improved evaluations of fi n t(n, <&, $*, T). In fact, we will 
see in Sec. V that such improvements can be appreciable 
even when the interactions are not particularly strong. 
Furthermore, <f>-DFT allows treatment of systems with 
significant inhomogencities, which are not describable by 
the simple two-fluid equations. 



Before closing this section, it is appropriate to state 
explicitly the differences in treatment which obtain for 
a DFT of bosons with only a single density. At zero 
temperature, one has nth = or n = $*$, a single 
density treatment would have fi n t{n) = gn 2 /2, and the 
Gross-Pitaevaskii equation [Eq. (|5"!?j) without nth] arises 
as the ground-state solution of the Schrodinger equa- 
tion, and need not be derived by shifting the quan- 
tum operator as in Eq. (I26|) . It is thus seen that 
in this limit the present treatment does not differ sig- 
nificantly from the single-density DFT treatment sug- 
gested by Nunes [37j]- Substantial differences do arise 
at finite temperatures, where a single-density treat- 
ment with a first-order local approximation would have 



fit 



(.9/2) ( 



2n 2 



in - C(3/2)A, 



and the corre- 



sponding effective potential, i> ff = v cxt -\-g (n + nth) with 
n th = C(3/2)A^ 3 , would still give rise to an equation for 
the ground state which is essentially the correct Gross- 
Pitaevskii equation, but the excited states would "feel the 
wrong potential" . In the limit of weak inhomogencities, 
the situation can be remedied. The Thomas-Fermi ap- 
proximation holds, with v n i — n = at points with a 
condensate, i.e., with n > £(3/2)A^ 3 . The correspond- 
ing free energy function is / (u) with u = A^n, and is to 
be continued to large densities, u > C(3/2). According 
to the rules for Legendre transforms, Eq. (|4T)) it is sim- 
ply linear in this regime. The solution of the Kohn-Sham 
system of equations for points with u n i(r) = n would then 
seem to be ambiguous, as there is a range of densities for 
a single value of the effective potential, but the condition 
i'ext(r) + Vi n t(n(r)) = fj, may be used to determine the 
density n(r) instead. If the interaction energy /; n t(n) 
correctly accounts for the difference between the refer- 
ence system and the interacting system, then the correct 
results for the free energy and the density distribution are 
guaranteed to obtain. It is only in the presence of signifi- 
cant inhomogencities that the weakness of this approach 
(i.e., the effect of its having essentially the same potential 
in the Schrodinger and the Gross-Pitaevskii equation) 
will show up. 



IV. DENSITY FUNCTIONAL THEORY FOR 
BOSONS WITH ANOMALOUS TERMS — A-DFT 

In this section, the thermodynamic approach will be 
used to develop another version of DFT for bosonic sys- 
tems, which results from adding a term of the form 

— J dr + £* tpi/j to the Hamiltonian. Here £(r) 

is a second fictitious potential — an anomalous potential 

— which is to be set equal to zero in the fully interacting 
system, £ cxt = 0. The fact that it does not vanish in 
the Kohn-Sham reference system, £ n j ^ 0, will result in 
a level of treatment generalizing that of Bogoliubov. In 
order to assist the reader, the partitioning into subsec- 
tions here is precisely parallel to that of the above section 
presenting $-DFT. 
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A. The grand potential and the free energy 

The thermodynamic treatment of the enlarged Hamil 
tonian follows the same steps as above, with Eq. (fTB] 
defining the grand potential, which acquires a [£, £*] de- 
pendence. The corresponding derivative is 



A(r) 



SQ 



(^(r)V(r)) , 



(53) 



where £(r) and £*(r) as well as rj(r) and 77* (r) play the 
role of B(r). The Legendre transform leading to the 
Hohcnberg-Kohn free energy is 

F H k(K$,$*,A,A*],T,A) = 

-/ rfr{ (,-,)n-^-^-^-rA }l 



and we have the additional relation 
SF HK 



SA 



(54) 



(55) 



B. The Bogoliubov Kohn— Sham System 



and should thus not be referred to as a thermal compo- 
nent. Note the complete cancellation of terms of type 
/ dr £* i < J' 2 in the effective Hamiltonian - the contribu- 
tions from $*(...)$ and £*;$ 2 are equal and opposite, 
due to Eq. (f57|) . In contrast, the terms of type / drrj ni Q* 
only partially cancel, leading to the \ prefactor in H con . 

The Hamiltonian H nc is quadratic in the field opera- 
tors, but does not conserve particle number. This form 
of Hamiltonian is diagonalizcd by the Bogoliubov trans- 
formation [53j]. Following Fetter's notation 54 j, the field 
operators may be written as 



(r) = Y^ u A v )lj - v * ( r )7j 



(60) 



and its Hermitian conjugate, where the primed sum runs 
only over positive energy solutions, £j > 0. The 7T and 
7j are bosonic creation and annihilation operators for the 
Bogoliubov excitations of the system. The generalized 
Schrodingcr equation for the (uj,Vj) wave functions is 
given by 



ft 2 V 2 
2m 

h 2 V 2 
2m 



+ v ni - fj, 
+ v ni - (1 



ij - ZZniVj = £j1 



(61) 



The Kohn-Sham reference system is described here by 
the nonintcracting Hamiltonian 



Hni = Jdrft (- 



h 2 V- 

2m 



dr 



{vni^ + Vni^ + ^^+CnM} , (56) 



where the nonintcracting effective potential i> n i and aux- 
iliary fields ?7 n i, ^ n i, are again to be defined by Eq. pT|) 
and determined by the interactions. Shifting the field op- 
erator by a scalar as in Eq. (j2"5"j) . $(r) = 3>(r) + 4>(r), and 
requiring all terms linear in the operators (ft (and <f) to 
cancel from the Hamiltonian, yields in this case 



ft 2 V 5 



2 m 



Vni - M * - Vm - 2£„i$* = . (57) 



The Kohn-Sham Hamiltonian becomes 

H n i = H nc + H con , (58) 

with 

H nc = [dr {^(-^- + v ni - 1*)$ - tniftft - CM , 

(59) 

and H con = -| / dr(r] eS $* + r}*^) as before, Eq. (1271) . 
The subscript nc represents the non-condensed part of 
the boson system, which persists to zero temperature, 



This two-component system of equations, known as 
the Bogoliubov-de Gennes equation [HH, is of the type 
H (") = £a z (") where % is a Hermitian matrix differen- 
tial operator, with inner product J dr (u* v*)a z (" J ) = 



fdri 



*i u 3 



1 

-1 



v*Vj) involving the Pauli matrix a z = 
The normalization is J dr (\iij\ 2 — \vj\ 2 ) = 1, 



and the orthogonality conditions are J dr (? 



i u 3 



for i ^ j and / dr (u*v* - v*u*) = for all i and j with 
Si,Sj >0@. 

Substituting Eq. (f6"0"|) into Eq. (j5"9")) . the Hamiltonian 
becomes 



H nc - ^2 £ J 



fj7j ~ / dr\ Vj \ 



(62) 



which is now in the form of a simple harmonic oscillator 
for each excitation mode j. With Eq. (|6"2")) . expectation 
values of different quantities at a temperature T can be 
evaluated, using either Eq. @DH with (7,-) = ("ft) = and 

(%"fj) = 6ij(exp(£j/T) — etc., or by explicitly cal- 

culating the partition function and the grand potential, 
and taking its derivatives. One finds that (V>(r)) = 3>(r), 
the density is n(r) = n nc (r) + $*$, with 



:(r) = 



\exp(£j/kftT) — 1 



(63) 
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and the anomalous density is A(r) = A nc (r) + <I> 2 , with 



Thomas-Fermi Approximation for A-DFT 



A„c(r) 



cxp(£j/kBT) — 1 



The grand potential of the Kohn— Sham system is fi r 

fine + ficon, With 



(64) 



finc(Ki - H, £ni, Ci] J T ) 



J2 (fcBTln 



1 — cxp 



Jj_ 

knT 



£j J dr\ Vj 



(65) 



and ficonbm-M^ni.C^ni.Ci] = ~1 I dr (Vni^* + 

as in Eq. . The Hohcnberg-Kohn free energy, accord- 
ing to Eqs. gjU and (J57J), is F ni = F nc + F con , with 

F nc ([n-<f>*$,A- $ 2 , A* - $* 2 ] ,T) = 
^ ffc B Tln(l - exp(-f J /fc B r)) - Ej J dr \ Vj A 

- y dr {?i nc (vni - m) ~ A nc £ ni - A nc £ ni } , (66) 



and F con [$,$*] = Jdr$* 



2m 



$ asinEq. The 



RHS of Eq. ([55]) is evaluated, as before, with the poten- 
tials «ni, Cni and £* ff which reproduce the non-condensate 
parts of the density and the anomalous density, through 
Eqs. ([6i"T) . (|63| and (|64|) . Notice that the non-condensate 
contribution to each of the thermodynamic quantities fi, 
n, A and -F can be further divided into a temperature- 
dependent thermal part and an athermal part, e.g., 
fi nc = fith + fiath in Eq. (|65[) . The temperature depen- 
dence yields an exponential convergence of the thermal 
parts, and only the athermal parts depend on the cutoff 
k c substantially (the above-mentioned rule for the evalu- 
ation of a Legendre transforms of a sum of two functions 
applies for F nc = F^ + -Path, with requirements such as 
SF t h/Sn th = 6F at h./5n a th and n th + n ath = n nc implied). 

In summary, the Kohn-Sham equations of A-DFT are 
Eq. (fS"T| for the condensate amplitude and Eq. (|6"Tj) for 
the non-condensate eigenstates and eigenvalues, together 
with the corresponding expressions for the density and 
the anomalous density, Eqs. (f6"5)l and (|6"4"|) . and together 
with the self-consistent determination of the effective po- 
tentials through Eqs. (|11[) and <jTUJ) . The interaction con- 
tribution to these potentials will be made explicit below. 
Once this system of equations has been solved, one may 
use the results to obtain the grand potential for the in- 
teracting system, which evaluates to 

fi = fine + J dr { /int - nv int + A£*„ t + A*£ in t + 
(7^** +tf nt *)/2} , (67) 
in full analogy with Eq. (gTJ) of $-DFT. 



For applications involving a large number of bosons, a 
Thomas-Fermi type of approximation can be formulated 
also in the presence of the anomalous potential. It is 
convenient to refer to a momentum variable p = Kk, with 
the density of states in the single-particle phase space 
taken as (2irh)~ 3 , as above. The corresponding "local" 
wave function, (itk, ^k) exp(zk-r), consists of plane waves 
with a "bare" energy of £k(r) = (h 2 k 2 /2m) + v n i(r) — fi 
(including the position dependence due to the effective 
potential). Eq. (|6Tj) then takes the form, 



1% _2Cni r =4. 2 r k 



(68) 



with the appropriate continuum normalization |uk| 2 
|t>k| 2 = 1- Solving this eigensystem of equations gives 



Uk = cosh 0k 7 
u k =(£y|£ ni |)sinh0 k 

£ k = £k/ cosh2^k = 



(69) 



[11 1 2 ) 



where tanh26*k = 2|£ n i|/£ k - 

When the non-condensate parts of the thermodynamic 
quantities are expressed in terms of the solutions of 
the Bogoliubov-de Gennes equation, they naturally have 
thermal and athermal parts, as in Eq. (|65[) , and in Eqs. 
([63]) and (|64|) . It will be convenient here to introduce di- 
mensionless functions / tn and / at h for the corresponding 
contributions to the grand potential, within the Thomas- 
Fermi approximation: 



fi„ c ~fc B TA^ / dr/ t „ 



EcK / dr/, 



/ath 



ksT ' ks,T J 



(70) 



where E c = h 2 k 2 /2m is the cutoff energy, |£ n i| is used as 
shorthand for */ £cff £ ni > an< ^ the functions are defined as 

f th (x, y) = —J q 2 dq In (l - cxp(- {q 2 + x) 2 - 4y 2 )) 

(71) 

and 



hih{x,y) = -j= J q 2 dq( K ^f{(i 2 



- V - cf 



(72) 

[the notation x, y, x and y will be used below as short- 
hand for the corresponding combinations in Eq. (|70j) ]. 
The grand potential of the Kohn-Sham system is thus 
the sum of three terms, a condensate part, a thermal 
part, and an athermal part. As the integral of the ather- 
mal part is divergent, we have used the cutoff scale to 
express it in dimensionlcss terms (the cutoff dependence 
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of the thermal part is exponentially small, and has been 
ignored). The cutoff energy is large, and it will be ap- 
propriate to expand / atn for small values of its variables 
— see the appendix. In contrast, the temperature may 
be small, and so the whole range of fth will be relevant, 
except in specific cases such as at T = 0. 

The square root in the integrand of Eqs. (|7Tj) and 
(|72[) is the normalized energy £k, and for x = 2y or 
"ni — A* = 2 |£ n i| it has a linear dependence at small 
wavcnumbers. This represents the phonon branch of the 
excitation spectrum of the superfluid. The present de- 
scription allows also for situations with x > 2y, which 
possess a gap in the spectrum at k = 0. As will be dis- 
cussed further below, this gap is not physical. 

The Thomas-Fermi results for the density and the 
anomalous density are: 



n nc (r) ~ X T 3 u th (x,y) + klu ath {x,y) , 

Anc(r) - 54 / (X^, 3 w t h(x, y) + klw ath {x, y) ) 

y ^ni 

with the following notation for the derivatives: 



(73) 



dfth(x,y) 
dx 
dfth{x,y) 
dy 



u a th(x,y) 



df at h(x,y) 



w at h{x,y) = - 



dx 

df a th{x,y) 



dy 



(74) 



The entropy has a contribution only from the thermal 
part, and is given by 



dr ( 5 

'- k B I jp- \ - -fth (x, y) + xu th (x, y) - yw th (x, y) 



The Hohenbcrg-Kohn free energy becomes 
F nc (K c ,A nc! A* c ],T)~ 

dr ^-^-/th (Uthi Wth) + Sckcfath (Math, Wath. 



(75) 



(76) 



where the Legendre transforms are defined as [again, the 
minimization may be replaced by the requirements of 
Eq. flU 



fth (u thl w t h) 



max 



(fth (x, y) - xuth + ywth 



/ath (Math, w a th) = max (/ ath (x,y) — xu at h + yw at h ) 

x,y 



(77) 



and are used with the conditions 



A^Mth + k 3 u a th = n nc , 2 A| wth + 2k 3 w a th = |A nc | 



k B T 



dfi 



Hi 



du 



£, 



Of 



at li 



th 



du 



, rrdfth c Of ath 
KBJ ~ = Or. 



ath 



dw 



th 



' dw, 



ath 



(78) 



Further details regarding / tn and / a th are given in the 
appendix, and for f t h also in Sec. V. 

As for the <3?-DFT of the previous section, a local, 
Thomas-Fermi approach can be applied to the conden- 
sate amplitude $ as well, amounting to neglecting the 
derivative term in Eq. (|57j). resulting in a completely lo- 
cal description. 



C. Interaction Effects 

A local description of the interaction effects (LDA) 
for the present application of DFT requires knowledge 
of /int(", A, A*, T) — the interaction energy den- 

sity for a uniform system of density n, condensate ampli- 
tude $, and anomalous density A, at temperature T. Its 
derivatives will determine the potentials of the reference 
system, which include 

dfint 



&i(r) = £ext(r) + £int(r) = _ 



OA* 



(79) 



in addition to Eq. (|40|) for the interaction contribution 
to the effective potential v n i and fictitious field r] n i. 

As in the previous section, one may use the weakness of 
the interactions in order to derive a simple approximation 
for /i n t, by evaluating the interaction energy to leading 
order in the interaction strength g. 



f 



int 



■4|$| 



(80) 



*2 



2,tlt\ 



f,txl\ 



due to the mixing of creation and annihilation opera- 
tors in Eq. (|60|) . since expectation values of the type 
{4>4>) n ° longer vanish. Similarly, Wick's theorem be- 
comes ($<^ 3 4 ) = ($$)(cfo0 4 ) + (4>\4> 4 )(4>l$3,) + 

{4>[h){4>U^)- Using both n = |$| 2 + (fifo, and A = 
$2 _ 



gives 

/int 



21*1 



(81) 



This result can be interpreted as a triple counting of 
the condensate term, compensated by a double sub- 
traction — the I "J 1 4 contribution appears in the direct 
((^i^4)('*/>2'03)), the exchange (($$3) {4>\^i)) and the 
anomalous ((^i^l) (^3^4)) term, but physically should 
be accounted for only once. Using this approximation in 
Eqs. (|4T)f and (f79")) for the effective potentials gives 

Uni(r) ~ Vcxt(r) + 2gn(v) , 
ry ni (r)~2 5 |$(r)| 2 $(r) , 



(82) 



Mr) 



■§A W 



Substitution of this approximation into Eq. (|57j) gives 
the generalized Gross Pitaevskii equation 

H 2 V 2 \ 
— — + Uext - M + 2gn nc + g\<S>(r)\ 2 ) $+gA nc <P* = , 
2m J 

(83) 
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which contains an extra term involving A nc [cf. Eq. (|52[)] . 
These potentials are also to be used in the Bogoliubov-de 
Gennes equation, Eq. (fBTj) . 

This approach reproduces the equations of the 
Hartree-Fock-Bogoliubov method, which have been sys- 
tematically derived and studied for a long time (l5l 151 - 
l48l |. In that context, the method is intended to calcu- 
late not only the thermodynamic properties of the sys- 
tem, but also its excitation spectrum. It has been crit- 
icized for producing a spectrum with a gap at small 
wavenumbers in homogeneous systems, while it is known 
that the correct long-wavelen gth result involves a gap- 
less, linear phonon spectrum j44|. A related difficulty 
arises at short distances, where it is seen that A(r) = 
lim r '_ H .(?/>(r)^>(r')) depends on the long-wavenumbcr 
cutoff k c (for T = and weak interactions, one may 
use the lim it of / at h discussed in the appendix to find 
n nc /n ~ \J nay and A nc /n ~ k c ao, in agreement with 
the literature). As A nc involves a product of opera- 
tors at essentially the same point in space, such a cut- 
off dependence should be accepted, and indeed, the na- 
ture of the point interaction should be expected to gen- 
erate a relationship between this ultraviolet divergence 
and long-wavelength behavior. Approximations which 
are designed to produce a spectrum without a gap have 
been studied [3, |47[, and more recently, a pseudopo- 
tcntial which allows for a rigorous treatment overcoming 
these difficulties (at least at zero temperature) has been 
suggested [5l| . 

In the context of DFT, it may be argued that the spec- 
trum is irrelevant, as it represents a property of the ref- 
erence system which need not be shared with the inter- 
acting system. In principle, it may thus be claimed that 
A-DFT is rigorously exact, despite the presence of the 
gap. However, it is clear that the spectrum affects the 
thermodynamic properties, and that if the reference sys- 
tem has properties which differ significantly from those 
of the interacting system, it will be difficult to find work- 
able approximations for F{ nt . It may thus be desirable to 
adopt the advanced pseudopotential approach to DFT 
[5l| . This would require using a reference system with 
nonlocal effective potentials, i.e., introducing terms of 
the form £(r, r')ipi(r)ip'(r') and v(r, r')tp' (r)ijj(r') into 
the Hamiltonian, with £ ox t = u xt = but £ n i,fni 7^ in 
close analogy to the derivation above, and is beyond the 
scope of the present work. 



V. THE UNIFORM-SYSTEM, 
WEAK-INTERACTIONS LIMIT 



modynamic perturbation theory. It will have relevance to 
the local density approximation, which (as noted above) 
relics on knowledge of the properties of uniform systems. 

The Thomas-Fermi approach is accurate (not an ap- 
proximation) for homogeneous systems. Hence, the re- 
sults of the Thomas-Fermi subsections above are directly 
applicable. Similarly, the lowest order, linear term in the 
interaction strength g has been given above, for both $- 
DFT and A-DFT. The discussion here pertains to a finite 
temperature, not too near to either the BEC transition 
or to T = 0, and thus, the interactions are treated as 
perturbing a noninteracting Bose-condensed system. 



A. Two— fluid method 

Applying the $-DFT of Sec. Ill to a uniform system, 
one drops the gradient terms in the Gross-Pitaevskii 
equation (|52p. and solves it in conjunction with Eq. (|46p . 
in the g —> limit. The external potential is set to 0, and 
the x — > limit of the thermodynamic function / defined 
in Eq. (gSJ is pertinent (cf. Fig. [2): 



/(ar) = -C(5/2) + C(3/2)x- 



giving 



4^ 



x 3 / 2 + 0(.t 5 / 2 ) , (84) 



u = f'(x) = C(3/2) - 2^x 1 ' 2 + 0(.t 3 / 2 ) , (85) 

where x = (2gn — /1) fk^T from Eq. (|51[) and u = A^n t h- 
This relation between the effective potential and the den- 
sity may be inverted as 



x = ^! + o (fe 4 ) 

47T V ' 



(86) 



where the notation Su = £(3/2) — u has been used, and 
the Legendre transformed Helmholtz free energy is, in 
dimensionless form, 



~ fin 3 
f(u) = -((5/2) + — + 0(Su 5 ) 

1Z7T 



(87) 



The Hohenberg-Kohn free energy per unit volume of 
the uniform system, including the interaction terms to 
leading order in g from Eq. (|49[) , is thus: 

F H K([n, T) k B T - 



V 



^-/(A£(n -*•*)) 



! (2n' -(*•*)' 



Having introduced both $-DFT and A-DFT, a discus- 
sion of applications is in order. Essentially all applica- 
tions are beyond the scope of the present article, but one 
"application" which is particularly revealing will be pre- 
sented here: the weak-interactions limit of homogeneous 
systems. Strictly speaking, this is not an application of 
DFT, and should be viewed instead as an exercise in thcr- 



Although it is straightforward to apply the self- 
consistent Kohn-Sham equations to this system, which 
amounts here to using, e.g., Eq. (|55"j) . a more physi- 
cally transparent discussion will result from following a 
minimum-energy path, closer in spirit to the Hohenberg- 
Kohn approach. The physical requirement that the ex- 
ternal auxiliary field vanish, r] ext = 9Fhk/9 < I'* = 0, will 
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thus be imposed by minimizing the free energy of Eq. 
with respect to **, at a given overall density n. From 
Eq. (|87|) . this amounts to minimization of 



9 



Su 3 



67r k B T 



(\ 3 T n -(,{?,/ 2) + 5uf 



(89) 



with respect to Su = Af,(|*| 2 - ri) + C(3/2), with 
Su > 0, but small. At g = 0, one finds Su = or 
|*o| 2 = n ~ C(3/2)Ay 3 , whereas for small g one finds 



l*nl 



This gives 



1*1 



2A; 



l*o| 



a result whose accuracy should be questioned, as dis- 
cussed below. Note that the particularly soft Su 3 be- 
havior of / (u), corresponding to the nonanalytic cusp in 
/ (x) at x = 0, has resulted in a sensitivity to interac- 
tions which is displayed by the sharp y/g dependence of 
the condensate amplitude (in other words, the position of 
a shallow minimum is easily changed by a perturbation 
which is sloped in that region). It is straightforward but 
not particularly illuminating to obtain additional ther- 
modynamic results at this level of approximation, e.g., 
to find the relationship between the chemical potential fi 
and the overall density n. 



B. Bogoliubov method 

Application of the A-DFT of Sec. IV to the present 
problem requires expanding the thermodynamic func- 
tions defined in Eq. (|70| at small values of their argu- 
ments [see the appendix, Eqs. (|A4j) and (|A5jl ]: 



/th(z,!/) = -C(5/2) + C(3/2)a; 
2^ 



3 



((a + 2yf 2 + (x- 2yf 2 ) + 0(x 2 ) , (91) 



where the requirement 0<y<x/2is used to drop O (y 2 ) 
contributions. A comparison with Eq. (|84[) is interesting 
already at this stage: clearly, setting the anomalous po- 
tential to zero, y = 0, reproduces the *-DFT result, but 
taking the Bogoliubov spectrum with a vanishing gap 
corresponds here to setting y = x/2. and results in an 
extra factor of \J~2 in the x 3 ^ 2 term. The expansion of 
/ath begins with 0(x 2 ) terms, Eq. (|A9|) . and therefore 
the athermal component does not affect the results to 
leading order in g. 

The scaled density and anomalous density are 



and 



dx 



C(3/2) - ^ [yjx + 2y + ^x - 2y) (92) 



Wth 



th 



dx 



2yfii(y/x + 2y-y/x-2y} , (93) 



and inverting these relationships gives 



1 

47T 



Su 



1 



th 



y ~ — SuthWth 
Sir 



(94) 



where again Su t h = C(3/2) — uth- The result for the 
Hohcnberg-Kohn free energy is 

/th(uth,wth) = ~C(5/2) + ^-(^8u 3 h +^w 2 h Su t h)+0(5uf h ) 

(95) 

The overall free energy, including the interaction terms 
to leading order, Eq. (|8ip . but excluding the athermal 
contribution, is: 



(90) F HK (n,*,**,A,A*,T) 



V 

k B T\ T 3 f th (A|(n - |*| 2 ), X 3 T | A - $ 2 

^(2n 2 - 



(96) 



|A| 2 -2|$| 



In the present case, both the fictitious potentials rj and £ 
must vanish, i.e., the free energy is to be minimized with 
respect to both $ and A. In dimensionless terms, this 
corresponds to minimization of 



1 (I 



2tt V3 



-Su 



1 



-Su th w 2 h 



k B TX 3 



(A r n - C(3/2) + Su th 



2 ) 



-2(A|n-C(3/2) + 5 Uth ) 2 ) (97) 



with respect to both Su t h and w t h (in principle, an arbi- 
trary phase factor could be associated with the term, 
reflecting the relative phase between A nc and * 2 , but the 
sign used here is clearly optimal). The minimization re- 
quires 



Su 



1 



ih 



Wth 



(A|n-C(3/2) + fe th + ^) , 



SuthWth 



^( 4 „-c(3/ 2)+ ^-^) 



(98) 



To leading order in g we may use Su t h — and 
Wth — in the RHS, resulting in Su t h — Wth/2 ~ 



(^T n - C(3/2j), or 



1*1 



|*o| 2 + A 7 



2ng 



*o| 



(99) 



together with a similar contribution to A nc (at this or- 
der, we have x = 2y). The interaction-dependent cor- 
rection to the condensate fraction is here a factor of 
y/2 smaller than the result of the previous subsection, 
Eq. (|M| . This discrepancy will be discussed in the next 
subsection. The present result, which includes the effects 
of the anomalous density through w, is in accordance 
with the literature [13] (recall that g ~ 47rS 2 ao/m and 
A T = v / 2Trh 2 /mk B T). 
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C. Comparison of the two— fluid and Bogoliubov 
Reference System 

It is at first surprising that the two methods discussed 
above lead to two different results for the leading-order 
correction to the condensate fraction of a uniform BEC. 
After all, both methods are based on a straightforward 
expansion in the small parameter g, within a thermody- 
namic framework which is in principle exact. The only 
source of error can be the neglect of terms in the inter- 
action energy, /i n t, which are of higher order in g. The 
next paragraphs display these terms explicitly. 

Clearly, the level of accuracy used in the A-DFT de- 
scription above can be imported into $-DFT, by choosing 
the appropriate form for /j nt (n, 4>, <£*). In fact, it is obvi- 
ous that simply dropping the anomalous term in Eq. (|97p , 
setting w t h = 0, reduces it to Eq. ([89]) . However, the 
proper procedure is to minimize over iy tn , which corre- 
sponds to imposing the condition £ = 0, i.e., to requiring 
the vanishing of the anomalous potential rather than the 
u>th contribution to the anomalous density. From the sec- 
ond line in Eq. ([55]), this leads (to leading order in g) 
to 



Wth 



16-rrg \^n-((3/2) 



k B TX^ 5u th 
Introducing this into Eq. (|97l) gives 



(100) 



1 . 
7T 5u 



9 



th 



k B T\3 



(X 3 T n - C(3/2) + Sua,) 



16tt 



g \ 2 (\ 3 T n-( 3/2 y 



k B T\ A T 



5u 



th 



(101) 



where terms up to second order in g have been retained 
(one can check a posteriori that the neglected terms are 
indeed small). Minimization of this, to leading order, 
amounts to the requirement that 



8v? 



ih 



27rg 
fc B TA| 



(A|n-C(3/2)) 



(102) 



which displays how dropping terms of second order in 
g leads to a doubling of the result for Su^ h . It is thus 
clarified that a weakly-interacting Bose-Einstein con- 
densed system is situated near a singular point, asso- 
ciated physically with a "completely full" thermal cloud, 
u th = C(3/2). At this special point, terms which are 
of second order in the weak interaction parameter g are 
not relatively small, because they are divergent, with the 
small quantity <5uth appearing in the denominator. 



VI. SUMMARY AND OUTLOOK 

The thermodynamic approach (summarized in Sec. 
II) provides a general method for generating DFTs for 



bosonic systems in thermal equilibrium at finite temper- 
ature, and has been used to derive the equations of $- 
DFT (Sec. Ill) and of A-DFT (Sec. IV). The different 
DFTs use as references different types of Kohn-Sham sys- 
tems, which are subject to different fictitious potentials. 
The reference systems have quadratic, soluble Hamilto- 
nians, and the interaction effects are to be included via a 
local-density approximation. The latter must be based 
on knowledge of the free energy of interacting homoge- 
neous systems, which are subject to the fictitious poten- 
tials. For $-DFT, knowledge of the free energy of a ho- 
mogeneous interacting system as a function of fx, T and 
rj is required in order to supply /i n t (ji, T), and for 

A-DFT a £ field must also be allowed for. This type of 
information is generally available only to leading order in 
the interaction parameter g. In this limit, our results for 
<I>-DFT generalize the two-fluid approach [4l[ to inhomo- 
geneous systems, and the results for A-DFT reproduce 
the Hartree-Fock-Bogoliubov model. As has occurred 
for electronic systems, we anticipate that the necessary 
results beyond the leading order will be generated us- 
ing quantum Monte Carlo techniques. Such techniques 
have been developed for Bose-condensed systems [561 ]. 
but as the fictitious potentials rj and £ break particle- 
number conservation, different variants of the techniques 
may need to be developed to meet this goal. 

It is of interest to compare $-DFT and A-DFT to the 
attempt to apply DFT to Bose-condensed systems made 
in Ref. [36[, using the standard approach to DFT rather 
than the thermodynamic one. In this reference, only two 
functional variables were used — n and $, as in $-DFT 
- but the reference or Kohn-Sham system chosen em- 
ployed a Bogoliubov type treatment, similar to that used 
in the A-DFT above. Correspondingly, the Hamiltonian 
of the reference system depends not only on the poten- 
tials used, but also on the (anomalous) density, which 
is to be calculated self-consistently. This represents a 
difficulty which the present thermodynamic derivation 
avoids. 

In comparing the two DFT methods, it was found that 
in the limit of homogeneous, weakly-interacting systems 
at finite temperatures, $-DFT (or the two-fluid method) 
does not correctly reproduce the leading-order correc- 
tion to the condensate density, and that this flaw can be 
corrected by including higher-order terms in /i nt . This 
is due to the extreme sensitivity of the corresponding 
energy-minimization problem: (a) terms of order g in 
the energy cause a shift in the minimizing value of $ 
which is of order y fg 1 and (b) the second order term in 
/int is divergent, having an O (y/g) denominator, thus 
contributing to O (g 3 / 2 ) in the energy and to a signifi- 
cant change in the O (^/g) correction to The higher- 
order terms in /; nt for $-DFT were in this case obtained 
from a first-order A-DFT calculation, which amounts to 
re-expressing the Hartree-Fock-Bogoliubov model as a 
minimization problem (Sec. V). 

The comparison just mentioned demonstrates the ad- 
vantages of A-DFT - it uses a reference system which is 
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much closer in its behavior to the fully interacting sys- 
tem, and therefore the approximation introduced is much 
less significant. It is reasonable to expect that this ad- 
vantage will be significant for inhomogeneous systems as 
well. Thus, obtaining the requisite data for homogeneous 
interacting systems as a function of both r\ and £ is called 
for. As an interim step, applications of $-DFT, for which 
strong-interaction corrections to /i n t could more easily be 
acquired, should also be considered. 
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Appendix A: Integrals for the Thomas— Fermi 
approximation of A-DFT 

In Sec. IV, Eq. (|70p . integrals corresponding to the 
thermodynamic functions fth and / a th were introduced, 
and in Sec. V, the need to evaluate these integrals in 
the limit corresponding to weak interactions arose. The 
details of the evaluations are presented here. 

For the thermal contribution, integration of Eq. ([7Tj) 
by parts gives 



fth (x, V) = 



sj(q 2 + x) 2 - 4y 2 



3^ Jo exp(i/(g 2 + x) 2 - 4y 2 ) - 1 
q 4 (q 2 + x)dq 



(q 2 + x) 2 - Ay 2 
The second factor is separated as 

gV + *> = q 2 - x+ 
(q 2 +x) 2 -4y 2 1 

If {x- 2y) 2 {x + 2y) 2 



(Al) 



2 \q 2 + x - 2y q 2 + x + 2y 



(A2) 



where the term q 2 — x diverges with q, and the remain- 
ing terms converge rapidly. This may be used to write 
fth = h + h , with the integrals 1\ and I2 involving the 
divergent and convergent terms, respectively. 

At small x and y, the square root term can be expanded 

as 



y/(q 2 + x) 2 - 4y 2 = q 2 + x + 0(x 2 ) 



(A3) 



where terms of order y 2 are included in 0(x 2 ) because of 
the limitation x > 2y > 0. One finds 



h 



-8 



q 2 +x 



J exp(g 2 + x) — 1 
-C(5/2) + C(3/2)a; + 0(x 2 ) 



(q 2 - x)dq + 0(x 2 ) 



(A4) 



(note cancellation of terms at order x 3 / 2 ) and 



h = 



-4 



{x~2yf 



{x + 2yf 



J a \q 2 + x — 2y q 2 + x + 2y 
^(( 2 ;-2 2/ ) 3 / 2 + (x + 2,) 3 / 2 ) +0 ( a; 2 ) 



dq + 0{x 2 



(A5) 



where in the last equation (q 2 + x)/(c q +x — 1) was ap- 
proximated by unity because the remaining factors are 
already of relatively high order in x, and are small when 
q 2 ^> x. These results are used in Sec. V, Eq. (j§ll . 

For the athermal contribution, one may rewrite Eq. 
(EU) as 



/ath (x,y) 
2 * 



--f+ 



y/(q 2 + x) 2 - Ay 2 -q 2 -xjq 2 + 2y 2 dq , 

(A6) 



where the integral can be continued to infinity without 
divergence, as can be seen by expanding the square root 
as 



^(q 2 +x) 2 -Ay 2 =q 2 + x -%+0 (x 3 ) 



Rcscaling q by \fx, and defining 
h (a) -- 



(A7) 



(q 2 + l) z - a - q z - 1 \ q z + 



dq , 
(A8) 

where a = Ay 2 fx 2 is a variable in the range < a < 1, 
gives 



/ath (x, y) 



4 25 5 / 2 

JLf + ^ h{a)+0 (x 3 ) . (A9) 



The integral h (a) vanishes at a — 0, and increases with 
a to a value of 8v^2/15 at a = 1. This value may be 
derived as 



VY + 2q 2 - q 2 - 1 1 q- 



1 



dq = 



Q 5 Q 3 Q r r- < x du , A , 
lim -^--^- + |+ / > /S( ti _2)_ , (A10) 
Q^oo 5 6 2 J 2 2 

where u = q 2 + 2 and the result obtains from the lower 
limit of the last integral. The derivative of h is given by 



h'(a) 



1 



<r 



dq 



(All) 



and decreases from h'(0) = n/4 to h'(l) = l/y2. As 
these values are within 10% of each other, a plot of h is 
very nearly a straight line, and is not included. It is easily 
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seen that h' (a) = 7r/4 — (ir/64) a + O (a 2 ) for small a, 
whereas h! (a) = l/y/2 - (l/16y/2) (1 - a) log (1 - a) + 
0(1 — a) near a = 1, where the logarithm arises from 
the 1/q behavior of the integrand for (l-a)<i| 2 Cl. 
This motivates the approximation 



h(a) ~ c\a + C20i 2 + cy,a i 



where the coefficients c\ 



■ C4Q; 4 + c (1 — a) log(l — a) , 



(A12) 
0.807 50, 



c 2 = -lis ( n + 3 A> - - 5 - 768 9 x 10 ~ 2 , c 3 = fjV^ 



7T ~ 2.775 x 10~ 3 , c 4 



1.666 3 x 10" 



T58) 7F 



'11 
. 10 



128 



V2: 



2.209 7 x 10~ 2 arc chosen 



so as to reproduce the calculated properties. Numerical 
integration of Eq. (|A8[) (evaluation of the integrand at 
large q requires some care) shows that the maximum er- 
ror in this approximation is less than 5 x 10 ~ 5 . As the 
cutoff scale is large, x and y typically are small, and such 
attention to the small-x-and-y limit is appropriate. 
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